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Time and space-resolved simulations are performed for the density profile and entanglement- 
entropy of trapped ultracold fermions, following near-adiabatic, moderately slow and instantaneous 
switching off of the trapping potential. The out-of-equilibrium dynamics of the Mott insulator is 
found to differ profoundly from that of the band insulator and the metallic phase, displaying the 
self-induced stability of an out-of-equilibrium quantum protectorate. 
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The experimental realization of ultracold gases of 
fermionic atoms in optical lattices is one of the major 
scientific breakthroughs of the past years [HUE)]. In ad- 
dition to the Pauli exclusion principle, the physics of such 
systems is governed by three distinct energy scales: the 
kinetic energy of the fermions, the potential energy due 
to the confining trap potential, and the fermion-fermion 
interaction energy. Various numerical [J, |5|, [y, [7| and ana- 
lytical [8( techniques predict that this interplay gives rise 
to a characteristic spatially varying density profile, dis- 
playing coexistence of metallic, Mott-insulator and band- 
insulator-like regions in different parts of the trap. Very 
recently, evidence for such phase-separated density pro- 
files in three-dimensional fermion gases has been obtained 
experimentally @, 0]. 

Most such investigations have been directed at the sta- 
tionary states of the many-fermion system, because of 
the similarity to the possible ground states of strongly- 
correlated many-electron systems in condensed-matter 
physics. Optical realizations of trapped fermions on a 
lattice, however, also allow one to study the time evo- 
lution of such systems, much more directly and eas- 
ily than in solid-state experiments, and in great detail 

[a ES EH, El El- The present paper makes predic- 
tions for such experiments. Specifically, we study the 
time evolution of the characteristic Mott insulator, band- 
insulator and metallic phases after rapid, moderate or 
near-adiabatic switching off of the trapping potential. 
This allows us to address two fundamental problems 
of many-body physics: l How does an out-of-equilibrium 
Mott insulator behave? and l how does the time evolution 
of a Mott insulator differ from that of a band insulator 
and of a metallic phase' 7 To answer these questions we 
present numerical predictions of the out-of-equilibrium 
density profiles corresponding to the distinct states real- 



ized in the trap. 

Before describing our methods and results, we recall 
that in a completely different part of physics a similar 
transition from static to time-dependent (TD) investi- 
gations is taking place: the study of entanglement in 
many-body systems. Entanglement in such systems is 
commonly studied from thepoint of view of its connec- 
tion to quantum criticality 14, ll5l[l6j or as a resource for 
quantum computation [171 . Ilo . Il9j ] . Such studies are fre- 
quently directed at properties of (pure or mixed) ground 
states or other stationary states of the system under in- 
vestigation. The time evolution of entanglement has re- 
ceived attention in the context of 
adiabatic quantum computation, but numerical studies 
frequently consider only the very particular time evolu- 
tion resulting from a quantum quench, and mostly focus 
on bosons or on pure spins. Here we present, in paral- 
lel with our results on the time evolution of the density 
profiles of an expanding cold fermion gas, also the cor- 
responding time evolution of the entanglement-entropy 
profile. 

We describe the trapped fermions by the Hamiltonian 
H(t) =Hq + V(t), with 

H = -t c \a c 3<? + U ^2 n iT n H + ^2v l h l , (1) 

(ij),a i i 

V(t) = -a(T)^2vihi. (2) 

i 

In Eq. (fT]), which describes a one-dimensional Hubbard 
model within an harmonic trap, (ij) denotes nearest 
neighbor sites and ni a = c\ a Ci a , with a =|, J,, is the 
local density operator expressed in terms of fermionic 
creation and annihilation operators. U is the on-site in- 
teraction and t the inter-site hopping (below taken to be 
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the unit of energy). The operator V(t) in Eq. ([2]), where 
hi = J2a^i^' controls the switching off of the parabolic 
potential = ki 2 , via the amplitude a(r), with the tem- 
poral boundary conditions a(0) = for the static trap, 
and a(r — > oo) — 1 for the completely switched off trap. 
Our explicit choice for a(r) at all times is 



o(r) - e(r) 
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The rate of the switch-off of the trap is thus determined 
by ro. We consider three cases: tq = + , tq = 300 and 
To = 460, representing sudden, intermediate and near 
adiabatic removals of the trap, respectively. 

The ground-state density profile of the Hamiltonian Ho 
is extracted using the Bethe-Ansatz (BA) local-density 
approximation (LDA) (25| to lattice-density-functional 
theory [26] . From the ground state we generate the time 
evolution within Time-Dependent Density-Functional 
Theory (TDDFT)^. We use the lattice version of 
TDDFT, introduced in [28| for the spin-independent 
case, which makes use of a spin-compensated, adiabatic 
BA-LDA (A-BALDA). Finally, the TD entropy profile is 
generated by means of the entropy LDA [2§| , again using 
an adiabatic description (i.e. the TD entropy functional 
depends locally in time and space on the TD density). 
These static and TD LDA-like approaches for the Hub- 
bard model are described in detail in the original ref- 
erences, where they have been tested and benchmarked 
against exact diagonalization, density-matrix renormal- 
ization and quantum Monte Carlo calculations. Such 
comparisons show that these LDA-like schemes have an 
accuracy of the order of a few percent for energies, par- 
ticle densities and entropies. Their very favorable com- 
putational cost allows one to study systems of hundreds 
of sites for any boundary condition, even in the absence 
of simplifying symmetries. 

Figure [1] shows a representative ground-state density 
profile and entanglement-entropy profile. Mott insulat- 
ing (M), band insulating (B) and vacuum (V) regions 
correspond to flat regions in the density profile and local 
minima in the entropy profile. These minima are sep- 
arated by metallic regions, which in one dimension dis- 
play Luttinger liquid (L) phenomenology. We now adopt 
this representative density profile as the initial state for 
the subsequent time evolution with the full Hamiltonian 
H(t) = H + V(t). The A-BALDA time evolution is 
performed using a predictor-corrector, split-operator al- 
gorithm, with the on-site effective potential computed 
in the mid-point approximation. Numerical convergence 
was further checked by halving A (30j . 

Panel a) of Fig. [2] illustrates the time-and-space re- 
solved density profile resulting from very slow (near adi- 
abatic) switching-off of the trap. At time r = the initial 
density profile is that of Fig. [TJ As expected for near adi- 
abatic switching, for very long times the density evolves 
towards the ground state of the unconfined system, which 



FIG. 1: (Color online) Density (thick line) and entropy (thin 
line) profiles for a chain with L — 100 sites and periodic 
boundary conditions, N = 60 fermions (with equal numbers 
of spin up and down ) with on-site interaction U/t = 8 and 
trapped by a static parabolic potential of curvature k/t = 
0.05. The five colored symbols indicate representative sites in 
the band insulating (B), Luttinger liquid (LI and L2), Mott 
insulating and near-vacuum (V) regions. In panels b) and c) 
of the other figures, the same colors refer to the same sites. 



in our case corresponds to a uniform distribution of 60 
fermions over 100 sites on a ring. 

Panel b) gives a more detailed view of the time evolu- 
tion of the five representative sites. The thin horizontal 
line indicates the uniform density distribution (n = 0.6) 
expected for the untrapped system. All sites ultimately 
attain this density, but in very different ways. Not unex- 
pectedly, the density at the metallic regions (labelled LI 
and L2) is rather fluid, and starts evolving towards the 
uniform density as soon as the trap is reduced. Similarly, 
the vacuum region (labelled V) gets filled up almost im- 
mediately. Interestingly, at intermediate times the vac- 
uum receives more fermions than would correspond to 
the uniform state. This overshooting, occurring roughly 
between r = 100 and r = 400, is characteristic of a 
nonequilibrium process. 

The Mott regions (labelled M), on the other hand, 
maintain their density fixed at that at r = 0, and re- 
main unchanged for long time before they finally start 
melting around r = 350. This is what one might intu- 
itively expect: while the metallic phase behaves more like 
a liquid, the insulating phase behaves more like a solid. 
However, the full picture is more complex than that: once 
the density at the originally metallic sites reaches 1 from 
above (LI) or below (L2), it, too, develops the character- 
istic Mott behavior and resists melting for an extended 
period of time. The corresponding transient flat regions 
are clearly visible in the curves labelled LI and L2 in 
panel b). Interestingly, the effect is more pronounced for 
approach from below than from above. 

The band insulator regions might have been expected 
to behave similar to the Mott insulator regions, but this 
expectation is not borne out by the data. Rather, the 
band insulator is as fluid as the metal, and starts evolving 
towards the uniform state as soon as the trap is beginning 
to be reduced. This points at a fundamental difference 
between band and Mott insulators: stabilization of the 
former is a consequence of the potential the fermions are 
exposed too (as it is in the case of the common band insu- 
lator in crystalline solids), while stabilization of the Mott 
insulator is due to particle-particle interactions, which 
are only indirectly influenced by changes of the trapping 
potential. Essentially, the Mott insulator self-stabilizes, 
while the band insulator requires an external potential 
to be stable. 
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FIG. 2: (Color online) Panel a): Time and space resolved density profile for adiabatic switching-off of the trapping potential. 
At time r = the density profile is that of Fig.[TJ while for later times the amplitude of the trap is slowly reduced, allowing the 
density profile to expand. Panel b): Cross sections of panel a), showing the time evolution of the density at the representative 
sites indicated in fig. [TJ Panel c): Time evolution of the entropy at the same sites. 



Panel c) displays the corresponding time evolution of 
the entanglement-entropy. Here we find that the two 
curves that displayed the most different behavior in the 
density (corresponding to vacuum and band insulating 
regions, respectively) display very similar behavior. The 
explanation of this reversal is that the entanglement- 
entropy is related to the degrees of freedom that are 
available for storing or recovering information, and this 
number is zero if a site is completely filled or completely 
empty. 

For longer simulation times the entropy also evolves 
towards that of the uniform system, but again this evo- 
lution is significantly retarded in the case of the Mott in- 
sulator, whose resistance to melting (increased disorder) 
is clearly seen to last until about t = 250. Unlike the 
total particle number, the total entropy is not conserved, 
and attains a maximum in the long-time limit. Since we 
are dealing with a nonequilibrium process, however, the 
approach to this maximum is not necessarily monotonic. 

Figure[3]presents similar results for moderate switching 
off of the trap. Here the approach to the final state is 
characterized by oscillations around the uniform state, 
which itself is not reached even for the longest simulation 
times available to us. However, the overall features of 
both the density and the entropy are still the same as for 
near-adiabatic switching. 

Finally, Fig. [4] shows results for instantaneous switch- 
ing. Here, the uniform state is never reached (hence, we 
display a much shorter time interval), and both the den- 
sity and the entropy profiles display the wild oscillations 
characteristic of a system driven forcefully out of equilib- 
rium. Even here, however, the Mott insulator resists the 
fluctuations more than the band insulator. 

In conclusion, we have found that the time evolution 
of the Mott insulator differs profoundly from that of the 
band insulator and the metallic phase. The difference 
reflects itself in self-induced stability of the Mott fea- 
tures even away from equilibrium, both in the density 
and the entanglement-entropy. This phenomenon is the 
dynamical equivalent of the quantum protectorate effect 
previously described in static situations [3 lj .and should 
reveal itself in time-of-flight experiments [2j, |7J . 

CV is supported by ETSF (INFRA-2007-211956). KC 
is supported by FAPESP and CNPq. 



[1] I. Bloch, J. Dalibard and W. Zwerger, Rev. Mod. Phys. 
80, 885 (2008). S. Giorgini, L. P. Pitaevskii and S. 



Stringari, Rev. Mod. Phys. 80, 1215 (2008). 
[2] R. Jordens, N. Strohmaier, K. Giinter, H. Moritz and T. 

Esslinger, Nature 455, 204 (2008). 
[3] M. Greiner and S. Foiling, Nature 453, 736 (2008). 
[4] M. Rigol, A. Muramatsu, G.G. Batrouni, and R.T. 

Scalettar, Phys. Rev. Lett. 91, 130403 (2003). M. Rigol 

and A. Muramatsu, Phys. Rev. A 69, 053612 (2004). 
[5] X.-J. Liu, P.D. Drummond, and H. Hu, Phys. Rev. Lett. 

94, 136406 (2005). 
[6] G. Xianlong, M. Polini, M. P. Tosi, V. L. Campo, Jr., K. 

Capelle and M. Rigol, Phys. Rev. B 73, 165120 (2006). 
[7] U. Schneider, L. Hackermiiller, S. Will, Th. Best, I. 

Bloch, T. A. Costi, R. W. Helmes, D. Rasch and A. 

Rosch, Science 322, 1520 (2008). 
[8] V.L. Campo, Jr. and K. Capelle, Phys. Rev. A 72, 061602 

(2005) . 

[9] W. Li, G. Xianlong, C. Kollath and M. Polini, Phys. Rev. 

B 78, 195109 (2008). 
[10] C. Kollath, U. Schollwock and W. Zwerger, Phys. Rev. 

Lett. 95, 176401 (2005). 
[11] C. H. Sc hunck, M. W. Zwierlein, A . Schirotzek, and W. 

Ketterle, |arXiv:cond-mat/0607298| 
[12] S. Iwai and H. Okamoto, J. Phys. Soc. Jpn. 75, 011007 

(2006) . 

[13] M. Rodriguez, S. R. Clark and D. Jaksch, J. Physics: 

Conference Series 99 012017 (2008). 
[14] A. Osterloh, L. Amico, G. Falci and R. Fazio, Nature 

416, 608 (2002). T. J. Osborne and M. A. Nielsen, Phys. 

Rev. A 66, 032110 (2002). 
[15] L.-A. Wu, M. S. Sarandy and D. A. Lidar, Phys. Rev. 

Lett. 93, 250404 (2004). L.-A. Wu, M. S. Sarandy, D. A. 

Lidar and L. J. Sham, Phys. Rev. A 74, 052335 (2006). 
[16] D. Larsson and H. Johannesson, Phys. Rev. Lett. 95, 

196406 (2005), ibid. 96, 169906(E) (2006). 
[17] M. A. Nielsen and I. L. Chuang, Quantum Computation 

and Quantum Information (Cambridge University Press, 

2000). 

[18] L. Amico, R. Fazio, A. Osterloh and V. Vedral, Rev. 

Mod. Phys. 80, 517 (2008). 
[19] V. Vedral, Rev. Mod. Phys. 74, 197 (2002). 
[20] A. M. Lauchli and C. Kollath, J. Stat. Mech. P05018 

(2008). 

[21] G. Sadiek, Z. Huang, O. Aldossary and S. Kais, Mol. 

Phys. 106, 1777 (2008). 
[22] L. Qiu, A. M. Wang and X. Q. Su, Opt. Comm. 281, 

4155 (2008). 

[23] M. Fagotti and P. Calabrese, Phys. Rev. A 78, 010306(R) 
(2008). 

[24] J. Fitzsimons and J. Twamley, Phys. Rev. A 72, 

050301 (R) (2005). 
[25] N. A. Lima, M. F. Silva, L. N. Oliveira and K. Capelle, 

Phys. Rev. Lett. 90, 146402 (2003). 
[26] K. Schonhammer, O. Gunnarsson, and R. M. Noack, 

Phys. Rev. B 52, 2504 (1995). 
[27] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 

(1984). 

[28] C. Verdozzi, Phys. Rev. Lett. 101, 166401 (2008). 



4 



FIG. 3: (Color online) As Fig. [2] but for moderate switching-off of the trap. 
FIG. 4: (Color online) As Fig. [2] but for instantaneous switching-off of the trap. 



[29] V. V. Franca and K. Capelle Phys. Rev. Lett. 100, 
070403 (2008). 

[30] The discontinuity in the effective single-particle potential 
321 . |33| ] makes the numerical time evolution challenging 
281 ] . as the discontinuity depends itself on the TD den- 
sity. To make the problem tractable, we (i) introduce a 
small smoothening of the discontinuity, and (ii) use a re- 
cursive time step in the time propagation, to ensure that 
in propagating between, say, t and t + A, no jump in v c 
was missed. After smoothing the discontinuity, the jump 
in v c broadens over a range Sn ~ 0.095, and its value 
is reduced by «s 10%. This modifies the shape of the 
Mott plateaus, which get slightly rounded (see Fig. [T}. 



Our conclusions about the Mott insulator resisting melt- 
ing longer than the band insulator are not affected. On 
the contrary, we expect that if the discontinuity is not 
smoothened, the Mott physics is even more pronounced 
and the plateaus will resist melting even longer. 
[31] R. B. Laughlin and D. Pines, P. Natl. Acad. Sci. USA 
97, 28 (2000). 

[32] D. Vieira, K. Capelle and C. A. Ullrich, Phy s. Chem. 

Chem. Phys. accepted (2009). larXiv:0902.1164l 
[33] N. A. Lima, L. N. Oliveira and K. Capelle, Europhys. 

Lett. 60, 601 (2002). 



This figure "ADIAsmall.jpg" is available in "jpg" format from: 



http://arxiv.org/ps/0905.1398vl 



This figure "INTERMsmall.jpg" is available in "jpg" format from: 



http://arxiv.org/ps/0905.1398vl 



This figure "STEPsmall.jpg" is available in "jpg" format from: 



http://arxiv.org/ps/0905.1398vl 



This figure "parabolacolor_small.jpg" is available in "jpg" format from: 



http://arxiv.org/ps/0905.1398vl 



